The melting curve of MgO from first principles simulations 
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First principles calculations based on density functional theory, both with the local density approx- 
imation (LDA) and with generalised gradient corrections (GGA), and the projector augmented 
wave method, have been used to simulate solid and liquid MgO in direct coexistence in the 
range of pressure < p < 135 GPa. The calculated LDA zero pressure melting temperature is 
T^ DA = 3110 ± 50 K, in good agreement with the experimental data. The calculated GGA zero 
pressure melting temperature T GGA — 2575 ± 100 K is significantly lower than the LDA one, but the 
difference between the GGA and the LDA melting curves is greatly reduced at high pressure. The 
LDA calculated zero pressure melting slope is dT /dp ~ 100 K/GPa, which is more than three times 
higher than the currently available experimental one ((Zerr and Boehler, Nature 371, 506 (1994)). 
As the pressure is increased, the melting curve deviates significantly from the experimental one. At 
the core mantle boundary pressure of 135 GPa MgO melts at T m — 8140 ± 150 K. 

PACS numbers: 64.10. +h 64.70.Dv 66.20.+d, 71.15.Pd 



MgO (periclase) is a very important material for a 
number of reasons. At ambient conditions it has the 
rock-salt structure of NaCl, and has the peculiarity of not 
showing any phase transition at least up to 227 GPa [1]. 
For this reason it is often used as a pressure medium 
in high pressure solid media devices. The periclase re- 
lated structure Mg(Fe)0-magnesiowustite is believed to 
be present in large quantities in the Earth's lower man- 
tle, therefore knowledge of the behaviour of MgO under- 
pressure is very important for our understanding of the 
lower mantle. In particular, the melting behaviour of 
periclase under pressure is important to put constraints 
on the solidus in the lower mantle: a low melting temper- 
ature of periclase may point to a low eutectic point, and 
in particular may support the suggestions of the presence 
of partial melt in the 'ultra-low velocity zone' [2]. 

The only available experiment for the melting be- 
haviour of MgO under pressure has been performed by 
Zerr and Boehler [3] (ZB), who measured the melting 
curve of MgO up to 32 GPa. They found a slope of 
the melting curve at zero pressure dT m /dp ~ 30 K/GPa, 
and a relatively low extrapolated melting temperature at 
lower mantle pressures; a result which may support the 
presence of partial melt at the bottom of the mantle. 

A number of theoretical calculations, based on empir- 
ical potentials, have attempted to determine the melting 
curve of periclase up to lower mantle pressures [4-8] . The 
results of these calculations are somewhat scattered, but 
they all predict significantly higher than experiments zero 
pressure melting slopes, and consequently much higher 
melting temperatures of periclase at lower mantle pres- 
sures. 

Here I report first principles calculations of the whole 
melting curve of MgO in the pressure range 0-135 GPa. 
Melting points have been calculated by performing di- 
rect first principles simulations of solid and liquid MgO 
in coexistence. Since a large number of atoms is needed 
to represent correctly solid and liquid in equilibrium, 



the method is extremely computationally intensive, how- 
ever, the feasibility of the coexistence approach within 
first principles calculations has been recently demon- 
strated [9-11]. In the constant volume constant inter- 
nal energy (NVE ensemble) approach to the coexistence 
method it has been shown that liquid and solid can co- 
exist for long times, provided V and E are appropriately 
chosen [9,12,13]. The average value of the pressure p 
and temperature T over the coexisting period then give 
a point on the melting curve. Size effects have been stud- 
ied quite extensively, and it was shown that correct re- 
sults, including in MgO, can be obtained in systems con- 
taining more than 500 atoms [6,9,14,15], although recent 
calculations on LiH have also been performed on system 
containing 432 atoms [10]. 

The zero pressure crystal structure of MgO is the same 
as the NaCl structure, and since no transitions have been 
experimentally observed up to at least 227 Gpa [1], I 
assumed that melting occurs from this structure. How- 
ever, I note that at high temperature this is not nec- 
essarily true, and melting may occour from a different 
crystal structure, as recently suggested by Aguado and 
Madden [16]. The existence of a more stable solid phase 
would increase the melting temperature. 

A possible source of error in the calculations may be 
due to the neglecting of the formation of defects in the 
solid. Although defects formation is not explicitly ex- 
cluded in the current approach, it is unlikely that the 
simulations are long enough to allow for significant ionic 
diffusion in the solid. A concentration c of defects would 
be responsible for a decrease of the solid Gibbs free en- 
ergy per atom g v ~ —k^Tc [17], and a consequent in- 
crease in the melting temperature 5T m ~ g v /s m , where 
s m is the entropy of melting. The concentration c of de- 
fects in MgO near the melting temperature is given by 
c = exp{ — E v /2kBT}, where E v is the formation energy 
of the defect. An estimate of c comes from the value 
of the formation energy of a Schottky defect, which in 
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MgO has a value between 4 and 7 eV [18]. Recent den- 
sity functional theory (DFT) [19,20] and quantum Monte 
Carlo [20] calculations point to a value close to 7 eV, but 
even using the lower value 4 eV, and T — 3000 K for the 
melting temperature, we obtain c~6x 10~ 4 , which re- 
sults in a completely negligible correction to the melting 
temperature. 

Present calculations have been performed using den- 
sity functional theory with various approximations for 
the exchange-correlation (XC) functional. I used the 
VASP code [21], with the implementation of an efficient 
extrapolation of the charge density [22], and the projec- 
tor augmented wave (PAW) method [23,24]. Single par- 
ticle orbitals have been expanded in plane-waves, with 
a plane-wave cutoff of 400 eV. The Mg potential has a 
core radius of 1.06 A, and the 3s 2 electrons in valence; 
the O potential has a core radius of 0.8 A and the 2s 2 2p 4 
electrons in valence. The structural properties of MgO 
in its rock-salt zero pressure crystal structure are com- 
pared with the experimental results in Table I. Pressure 
against volume curves are compared with experiments in 
the inset of Fig. 1. In Table I I also report the value 
of the transition pressure between the rock-salt and the 
CsCl structures, which is in the region of 500 GPa. I also 
tested a " small core" Mg potential with both 3s 2 and 2p 6 
electrons in valence. Results are reported in Table I. The 
zero temperature MgO crystal has a band gap of several 
eV's, but this gap is significantly reduced at high tem- 
perature, where electron excitations become important, 
therefore finite temperature calculations have been per- 
formed by minimising the Mermin functional [25]. In 
fact, the electronic entropy contribution to the free en- 
ergy of the system is non-negligible, and at zero pressure 
it lowers the free energy of the liquid with respect to 
the free energy of the solid by almost 0.1 eV. Given a 
zero pressure entropy of melting is ~ 2.2 /ce/atom (see 
Tab II), the electronic entropy is responsible for a lower- 
ing of the zero pressure melting temperature of ~ 500 K. 

The coexistence simulations have been performed us- 
ing the local density approximation (LDA) for the XC 
functional, and the Mg potential with only 2s 2 in va- 
lence, using the NVE ensemble for systems contain- 
ing 432 atoms (3x3x6 cubic supercell). The time 
step was 1 fs and the convergency threshold on the to- 
tal energy 2 x 10 -8 eV/atom. With these prescrip- 
tions the drift in the micro-canonical total energy was 
less than 0.5 K/ps. Simulations were carried out for 
up to 20 ps, and performed using the Y point only. 
Spot checks with Monkhorst-Pack [26] (2x2x1) and 
(2x2x2) k-point grids showed energies converged bet- 
ter than 0.1 meV/atom and pressure converged better 
than 1 MPa. A correction term of 3.7 (5.7) GPa due to 
the lack of convergency with respect to the plane wave 
cutoff has been added to the calculated pressures in the 
low (high) pressure regions. 

To prepare the system I have followed the same pro- 



cedure employed in Refs. [9,13]. A perfect crystal is 
initially thermalized to a guessed melting temperature 
T guess, then the simulation is stopped, half of the atoms 
are clamped and the other half are freely evolved at very 
high temperature until melting occurs, then the liquid is 
thermalized back at T guess . At this point the system is 
being freely evolved in the NVE ensemble. Note that 
in (V, E) space melting is represented by a "band" , and 
not by a line as in (p, T) space. In particular, for every 
fixed V, a whole piece of melting curve can be evaluated 
by varying E appropriately. The amount of total energy 
E given to the system can be tuned by assigning differ- 
ent initial values to the velocities [27]. Now, unless the 
value of E is in the right range for the volume V, the 
system will either completely melt or completely solid- 
ify. For MgO this happens very quickly, typically in less 
than 1 ps. Therefore, for each fixed V , a certain number 
of "trial and errors" steps are required in order to find 
the right value of E, for which solid and liquid coexist 
for long time. 

Simulations have been performed at a number of points 
in the pressure range 0-135 GPa. The points on the 
melting curve obtained from these simulations are dis- 
played in Fig 1, and compared with the experimental 
results of ZB, and with previously calculated melting 
curves [4-8]. The solid curve between points has been 
obtained by interpolating with third order polynomials, 
with the conditions for the coefficients given by the co- 
ordinates and the slopes at the points. The slope of the 
melting curve is obtained from the Clausius-Clapeyron 
relation dT/dp — v m /s m , where v m and s m are the vol- 
ume and the entropy change on melting respectively. For 
a chosen point on the melting curve (p m ,T m ), v m is cal- 
culated from independent simulations on solid and liquid 
MgO, with the respective volumes adjusted until the cal- 
culated pressures are equal (within ~ 0.5 GPa) to the 
chosen value p m . From these simulations it is also pos- 
sible to obtain s m , given by s m = (e m + p m v m )/T m , 
where e m = en q — e so i, with ey q , e so i the internal ener- 
gies in the liquid and the solid simulations respectively. 
These simulations have been performed on cells contain- 
ing 64 atoms, for over 40 ps, and spot checked with sim- 
ulations performed with 216 atoms, which showed essen- 
tially no difference in e m within a statistical error of ~ 10 
meV/atom. Pressures calculated with the 64-atom and 
216-atom cells differed by less than 0.5 GPa. Results for 
v m , Sm and the melting slope are reported in Table II. 

As a by-product of the simulations on the liquid state 
I have also obtained the shear viscosity of the liquid T], 
calculated using the Green-Kubo relations as described in 
Ref. [28], this is also reported in Table II. Note that the 
value of the shear viscosity on the melting curve increases 
only slightly as a function of pressure. 

In order to test the effect of the size of the simula- 
tion cell on the melting curve, I have also performed one 
simulation with 1024 atoms (4x4x8 cubic supercell) 
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at p ~ 47 GPa for 11 ps [29]. The melting point (p,T) 
extracted from this simulation is also reported in Fig. 1, 
and is essentially indistinguishable from that obtained 
with the 432-atom cell. 

Other systematic sources of errors are associated with 
the XC functional. In order to test how different ap- 
proximations for the XC behave, I have used the gener- 
alised gradient corrections (GGA)'s known as PW91 [30] 
and PBE [31], and evaluated GGA - LDA differences in 
the melting temperature at a number of different pres- 
sures. To do so, I have calculated the free energy differ- 
ences between LDA and the GGA functionals employed, 
for both solid and liquid, following the techniques de- 
scribed in Ref. [13]. At zero pressure the difference in 
free energy between LDA and PW91 is 0.1 eV/atom, and 
the difference between LDA and PBE is essentially the 
same. When combined with the entropy of melting of 
~ 2.2 fee /atom, this difference results in a lowering of 
the melting temperature of ~ 540 K. This result is in fair 
agreement with the findings of Tangey and Scandolo [8] , 
who also reported that at zero pressure the GGA melt- 
ing temperature is lower than the LDA by ~ 450 K. So, 
similarly to the case of Al [9,17,32], in MgO the GGA un- 
derestimates significantly the zero pressure melting tem- 
perature. However, this result should not be taken as 
a general trend of GGA versus LDA, as in Si, for ex- 
ample, the situation is reversed, with the LDA the zero 
pressure melting temperature being lower than the GGA 
one [33,34]. At higher pressures the free energy differ- 
ences between GGA and LDA are greatly reduced, and 
at 135 GPa the GGA melting temperature is lower than 
the LDA one by only ~ 100 K. 

Finally, I have tested the effect of the choice of the dis- 
tribution between valence and core electrons in the Mg 
potential, by evaluating the free energy differences be- 
tween the MgO system represented with the couple of 
Mg and O PAW potentials described above and the sys- 
tem in which the Mg PAW potential was chosen to have 
both 3s 2 and 2p 6 in valence. The free energy differences 
are extremely small, both at zero and at high pressure, 
and result in corrections to the melting temperature of 
26 K and 69 K at and 135 GPa respectively. These 
corrections have been included in the results reported in 
Fig. 1 

The present results for the melting curve of MgO, 
based on first principle simulations of direct coexistence 
of solid and liquid, support the "high" melting curve pre- 
viously indicated by a number of theoretical approaches 
based on empirical potentials [4-7], but are at variance 
with the experimental results of ZB. The agreement be- 
tween the LDA and GGA predictions for the high pres- 
sure melting curve of MgO suggests that the expected 
systematic error due to the XC functional employed may 
be small. However, this is not the case in the low pres- 
sure region, where the difference between LDA and GGA 
is of the order of ~ 18%. This relatively large discrep- 



ancy points towards the need of going beyond DFT with 
the current implementations of the XC functionals. We 
believe that one possible way forward will be the use of 
quantum Monte Carlo techniques [35]. 

The predicted high melting curve of MgO has impor- 
tant consequences for our understanding of the Earth's 
lower mantle and the history of the Earth's formation. 
The melting temperature of the mantle is constrained 
between the melting temperatures of the end members 
Mg(Fe)0-magnesiowustite and Mg(Fe)Si03-perovskite, 
and the eutectic point. The very high melting temper- 
atures of both MgO and MgSi0 3 (between 7000 and 
8500 K [36]) indicate that the eutectic temperature of 
the lower mantle is much higher than its current temper- 
ature (estimated to be between 2550 and 2750 K) and 
than the temperature at the top of the core ( ~ 4000 K 
). This would suggest that the presence of partial melt 
of Mg-bearing phase in the ultra-low velocity zone at the 
bottom of the Earth's mantle is unlikely. 
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FIG. 1. Melting curve of MgO obtained with present 
DFT-LDA coexistence simulations performed on 432-atom 
cells (black dots and heavy solid line), 1024-atom cell (black 
square), and present DFT-GGA results (triangles), compared 
with experiments (open diamonds) [3] and theoretically de- 
termined curves based on empirical potentials: Belonoshko 
and Dubrovinski (chained line) [6], Strachan et al. (dashed 
line) [7] , Cohen and Zong (dotted line) [5] , Vocadlo and Price 
(solid line) [4], Tangney and Scandolo (heavy dashed line) [8]. 
Inset: pressure as function of atomic volume for rock-salt 
MgO calculated with DFT-LDA (solid line) and DFT-GGA 
(dashed line), and compared with experiments (stars) [1]. 
Calculations do not include zero point motion. 
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do 



(A) 



Bo(GPa) 



P tr (Bl-B2) (GPa) 



> 227' 
503 
505 
491 



Experiments 
DFT-LDA (large core) 
DFT-LDA (small core) 
DFT-GGA (PW91) 



4.213" 4.211" 4.212 c 4.19 d 
4.151(4.180) 
4.165(4.194) 
4.234(4.263) 



160±2 a 160. T 156 e 164.6" 
180(170) 
177(167) 
158(148) 



TABLE I. Experimental and calculated lattice parameter a and bulk modulus B of MgO in the NaCl structure (Bl), and 
transition pressure P tr (Bl-B2) between the NaCl and the CsCl (B2) structures. Values in parenthesis include zero point motion 
and room temperature effects estimated from Ref. [37]. The calculations labelled "large core" have been performed with the 
Mg potential with only the 3s 2 electrons in valence, while those labelled "small core" with the Mg potential with both the 3s 2 
and the 2p 6 electrons in valence. 



Pm (GPa) 


T^ DA (K) 


T GGA (K) 


dT/dp (K/GPa) 


v m (A /atom) 


s m (fcB/atom) 


rj (mPa s) 


-0.4(2) 


3070(50) 


2533(100) 


102(5) 


3.08(5) 


2.19(10) 


3.0(4) 


17.0(2) 


4590(50) 


4405(100) 


62(3) 


1.44(5) 


1.69(4) 


3.8(3) 


47.0(2) 


6047(40) 


5887(90) 


33(2) 


0.73(3) 


1.59(3) 


4.5(4) 


135.6(2) 


8144(40) 


8044(80) 


16(1) 


0.34(2) 


1.51(2) 


5.0(4) 


TABLE II. 


Calculated melting properties 


of MgO: pressure p m , 


LDA and GGA melting 


; temperatures T^ DA 


and T GGA , 



slope of the melting curve dT/dp, volume and entropy change on melting v m and s m , and shear viscosity of the liquid r\ at the 
melting point. All quantities (except T GGA ) have been calculated using DFT-LDA. 



5 



